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Abstract 

We obtain an implicit equation for the correlation dimension which describes clustering of 
inertial particles in a complex flow onto a fractal measure. Our general equation involves 
a propagator of a nonlinear stochastic process in which the velocity gradient of the fluid 
appears as additive noise. When the long-time limit of the propagator is considered our 
equation reduces to an existing large-deviation formalism, from which it is difficult to 
extract concrete results. In the short-time limit, however, our equation reduces to a solv¬ 
ability condition on a partial differential equation. In the case where the inertial particles 
are much denser than the fluid, we show how this approach leads to a perturbative expan¬ 
sion of the correlation dimension, for which the coefficients can be obtained exactly and 
in principle to any order. We derive the perturbation series for the correlation dimension 
of inertial particles suspended in three-dimensional spatially smooth random flows with 
white-noise time correlations, obtaining the first 33 non-zero coefficients exactly. 
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1 Introduction 


In aerosols and other suspensions of microscopic bodies it may be satisfactory to neglect 
hydrodynamic interactions, and to assume that the particles move independently. It is 
known that small particles moving independently in an incompressible turbulent or com¬ 
plex flow may show a pronounced tendency to cluster. This occurs if the time scale for 
viscous damping, r p , is comparable to the smallest characteristic time scale for fluctuations 
in the flow, r. Maxey |l] proposed that these ‘inertial particles’ cluster because they are 
expelled from vortices by the centrifugal effect (if they are denser than the fluid in which 
they are suspended; bubbles are expected to congregate in vortices). Later, Sommerer and 
Ott [2] showed that, in common with other chaotic dynamical processes, the trajectories 
of particles advected on compressible surface flows approach a fractal attractor (Ott [3] 
gives a good introduction to the role of fractals in dynamical systems). Numerical ex¬ 
periments by Bee [5] confirmed that fractal clustering is observed for inertial particles in 
incompressible flows, just as in the compressible surface flows considered in [ 2 ]. 

This clustering is of fundamental importance to understanding the effect of turbulence 
on aerosols, because of its potential relevance to the coalescence of cloud droplets into rain 
[ 5 ], or of dust grains into planetary precursors [ 6 ]. 

The clustering process and its fractal dimension have been investigated numerically 
in many works: [7| and 0 report state-of-the-art contributions. The present theoretical 
understanding of this effect is reviewed in [9], 

The present paper is concerned with the analysis of the correlation dimension, D 2 , 
which is the most important dimension in physical applications, but which is still quite 
poorly understood. The importance of the correlation dimension arises from its direct 
relation to the two-point correlation function of particles (given as equation ([S]) below), 
which enters in theories for collision processes pm nu and light scattering }12l . 

Our approach gives an implicit equation for the correlation dimension, in terms of a 
propagator for a nonlinear stochastic process in which components of the velocity gradient 
tensor appear as additive noise. In the limit as the propagation time approaches zero, 
our equation becomes a solvability condition for a linear partial differential equation. 
We analyse this system using perturbation theory. We obtain a series expansion for the 
correlation dimension of the particle distribution in powers of a dimensionless parameter 
which measures the importance of inertial effects. 

We briefly review the state of the theoretical knowledge concerning the clustering of 
particles. Maxey’s original work [L| proposed that the particles (which we assume are 
much denser than the fluid) are expelled by centrifugal forces from vortices in the fluid, 
but that this effect can only be effective when the motion of the particles relative to the 
fluid is neither too lightly damped nor too heavily damped. The damping is characterised 
by a dimensionless number termed the Stokes number, defined by 

St = — 

7 r 

where r is a characteristic time scale of the fluid flow, and where 7 = l/r p is the rate 
constant for damping the motion of the particles relative to the fluid. He also showed that, 
when inertial effects are weak, the particle velocity may be approximated by an effective 
velocity field which has a compressible component. Simulations do show that particles 
have lower density in regions of high vorticity, see PI. 
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The particle distribution has clustering properties which are much more significant 
than the instantaneous negative correlation between density and vorticity. The particles 
approach a fractal measure. This can be characterised in a variety of ways, but the 
approach which is most easily understood and most fundamental to physical applications 
is to consider the number of particles M inside a ball of radius 5r centred on a randomly 
selected test particle. For sufficiently small values of dr, the average of this quantity has 
a power-law dependence upon 5r with exponent denoted by D 2 - 

(Af{6r)) ~ 5r° 2 . ( 2 ) 

Throughout this paper the expectation value of X is denoted by ( X }. The exponent D 2 
is termed the correlation dimension of the particle distribution [3] • The fractal dimension 
of particle clusters has been investigated numerically, and it has been confirmed that the 
fractal dimension in turbulent flows is significantly less than the space dimension only 
when the Stokes number is of order unity (see for example Fig. 2 in [L|j, and Fig. 1 in [8]). 

However, a theoretical analysis leading to quantitative results concerning the depen¬ 
dence of the dimension D 2 upon the Stokes number is lacking. There are a few works 
in which analytical results on the correlation dimension have been obtained. Most of the 
literature has discussed the relation between the Renyi dimensions and the statistics of 
the finite-time Lyapunov exponent: these relationships were established by Grassberger 
and Procaccia m (see also p76j ), and are reviewed in the book by Ott [3|. Usually the 
finite-time Lyapunov exponent can only be investigated numerically, but its statistics can 
be obtained for the Kraichnan model m in which a particle is advected in a velocity field 
with white-noise temporal correlations: Falkovich et a 1 [LB] discussed the calculation of the 
Renyi dimensions for the Kraichnan model. The first analytical studies on the correlation 
dimension for inertial particles was made by Bee et a 1 EE], who considered a velocity held 
which has white-noise temporal correlations. Their method yields the first two terms of 
the series expansion of D 2 , but it seems to be very difficult to extend to higher orders (and 
the second-order coefficient in m appears to be incorrect). In [20] we described a new 
method which related the correlation dimension to the solution of a partial differential 
equation. It was shown that the series expansion of the solution to this equation can be 
automated, so that coefficients of arbitrary order are obtained by repeated application of 
a system of annihilation and creation operators. In this way the coefficients in a series 
expansion of the correlation dimension of inertial particles in two-dimensional random 
flows was obtained. Gustavsson and Mehlig [21] used a different technique to compute 
the correlation dimension for a random-flow model in one dimension where the correla¬ 
tion dimension could be treated as a small parameter. In a series of papers Zaichik and 
Alipchenkov [22* i23i [24] developed an approach to calculating the clustering and collision 
rates of particles in a turbulent flow which combines empirical data on turbulence, with a 
stability analysis of the dispersion of particles. 

In this paper we describe a general principle (section [2j) for calculating the correlation 
dimension, based on the invariance of the distribution of small separations under dila¬ 
tions, corresponding to translations in logarithmic variables. In section [3] we show how 
this principle can be expressed in terms of a time-propagator. We show that a large-time 
expansion of the propagator gives a set of equations closely related to equations derived 
from a large-deviation principle - discussed in mmm- We also show how an approx¬ 
imate expression for D 2 can be recovered from the large deviation formalism, but it is 
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difficult to extend this because of the intractability of determining the entropy function of 
the large deviations of the Lyapunov exponent. A short-time expansion of the propagator, 
by contrast, yields a partial differential equation involving D 2 which is more amenable 
to analysis. This approach was previously outlined in [[20]. Here, in section [4] we apply 
the method to a white-noise random-flow model in three spatial dimensions, developing 
a perturbation theory for D 2 in section [5] Because the correlation time of the flow is 
r = 0 for our model flow, the Stokes number is not defined for our model. However, our 
perturbation parameter, e, plays a role which is analogous to St. The relation between e 
and St is discussed carefully in [25], where it is argued that e 2 oc St. The perturbation 
series is divergent and the methods used to extract finite results are discussed in Section 
El Section [7] contains our conclusions and discusses possible extensions of this work. 


2 The correlation dimension 


In Secs. [2] and El we define the correlation dimension and discuss several distinct but 
interconnected approaches to calculating it. What these approaches have in common is 
that they use a dynamical variable, Z\ (t), which is derived from the linearised equation of 
motion. The statistics of Z\ (t) are also closely related to the leading Lyapunov exponent. 
Several different probability density functions must be introduced. We denote the prob¬ 
ability density function (PDF) of a quantity X by a function px , so that the probability 
element for X to lie in the interval [X, X + dAl] is d P = px(X)dX. The expectation value 
of X is denoted by (X). 

The correlation dimension D 2 is defined in terms of the expected number ( J\f(dr )) of 
particles inside a ball of radius dr surrounding a test particle: 


D 2 = 


lim 
< 5 r —>0 


ln(AA(dr)) 

ln(dr) 


(3) 


so that 

{N{dr)) ~ dr° 2 (4) 

which is the volume element of a ball in D 2 dimensions. If D 2 = cL (where d is the 
dimensionality of space) there is no clustering. The probability density p(5r) for a particle 
to have another particle at small distance 5r is 


p(6r) 


d(A f{5r)) 
d<5r 


dr D2 1 . 


(5) 


Note that this quantity is the ‘two-point correlation function’ which plays an important 
role in physical kinetics mm and scattering theory [12] . 


2.1 Logarithmic separation dynamics 


It is not immediately clear why the limit in equation Q should exist. In this paper we 
show why it does, and how to extract information about D 2 by considering a quantity 
Zi(t) defined by 


dr- 

dr 


= Z 1 . 


( 6 ) 


4 




Here dr denotes the time derivative of 5r, and Z\ is the logarithmic derivative of dr. We 
also consider the variable 

Y(t) = ln<5r(i). (7) 

The two variables Y and Z\ are related by 

Y(t) = y(0)+ fdt'Z^t'). (8) 

Jo 

We will argue that, in the limit as Y(t) —>• —oo, the variable Z\ obeys an equation of 
motion which is independent of Y. This implies translational invariance in the statistics 
of Z\. Correspondingly the PDF pv(Y) of Y exhibits translational invariance: py(Y) and 
Py(Y — Yq) must be the same function, up to a normalisation factor, for any choice of the 
displacement Yq. Hence 

Py(Y) = C(Y 0 ) P y(Y-Y 0 ) (9) 

for some choice of C(Yq) . The solution of this equation is 

Py(Y) = Aex.p(aY) (10) 


for some constant a and normalisation A. This expression is valid only for Y —» —oo, 
so that we must require a > 0 to give a normalisable probability density. Consider 
the corresponding PDF of dr, denoted by p$ r (dr): the element of probability is dP = 
PSr(dr) ddr = py(T)dT = Adr a ~ 1 ddr, so that the distribution of dr corresponding to 

m ^ 

P5r(dr) = Adr^ 1 . ( 11 ) 

By comparison with ([5]) it follows that the exponent of the distribution of Y and the 
correlation dimension are equal: 

D 2 = a . ( 12 ) 

Thus we conclude that D 2 can be determined by studying the statistics of the logarithmic 
derivative Z\ = dr/dr. Specifically, if Z\(t) is a random variable with statistics that 
become independent of Y as Y —> — oo, then the distribution of Y is Py(Y) ~ exp(D 2 ^). 
So, to determine D 2 we need to study the equation of motion for Z\ (t) and how the 
statistics of Z\ determine the exponent a. 

Before going on to consider the equation of motion for Z\ 1 we remark that the vari¬ 
able Z\ (t) also gives information about the leading Lyapunov exponent A: provided the 
separations remain sufficiently small, we have 


A = lim 1/ln ( 

o t \ y<5r(0) J 


We can express this in terms of a limit of a finite-time Lyapunov exponent a(t): 

1 rt 
t Jo 

The leading Lyapunov exponent is therefore an expectation value of Z\(t): 


=K 1 " (£1) 


fdt' Z x (t') 

Jo 


A = lim a(t) = (Zi(t)) 

t—t OO 


(13) 


(14) 


(15) 
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2.2 Equation of motion for the logarithmic derivative 


We have shown that information about D 2 is contained in the dynamics of the logarithmic 
derivative of the separation, Z\(t). To proceed further we need an equation of motion for 
this quantity. The equations of motion for a small spherical body moving in a viscous 
fluid are discussed in [26, 27j- We consider the case where the density of the body is much 
higher than that of the surrounding fluid. In this limit the equations of motion for the 
particle position r(t) and velocity v(t) are: 

r = v, h = 7 [u(r(t),t) — v] . (16) 

An equation of motion for Z\ is derived from the linearised equations of motion describing 
a pair of particles with a separations dr and dv in position and velocity 


dr = 5v , dv = —7 5v + 7 E dr . (17) 

Here E is the matrix of flow-velocity gradients with elements E{j = dui/drj. From these 
equations we must obtain an equation of motion for Z\ = dr/dr, where dr = |<5r|. To illus¬ 
trate the approach in its simplest context, we show how this is done for a one-dimensional 
model, where x is the coordinate of the particle. In one dimension we have dr = |&c|, and 
simple manipulation of equations (fT71) gives 

Z\ = — 7 Z 1 — Z\ + 7 E(t) (18) 


where 


_. , du , , . . 

E (t) = • 


(19) 


In two or three dimensions, the variable Z\{t) is coupled to one or more additional vari¬ 
ables, but there are always a finite number of variables, Z\, Z 2 , ■ ■ ., which are coupled in 
a closed system of equations analogous to (|T8l) . 

The one-dimensional version of equation (11711 allows particles to exchange positions, 
that is dx passes through zero while dv remains finite. This corresponds to a ‘caustic’ 
singularity [28] where Y(t) goes to —00 and returns, while Z\(t) goes to —00 and returns 
from + 00 . This divergence of Z\ is a special feature of the one-dimensional version of the 
model and it is absent in higher dimensions. We should nevertheless consider its effect. 

The finite-time singularities give rise to a ‘tail’ of the distribution of Y. Consider the 
form of the distribution of Y resulting from a fold event in a one-dimensional system, 
where one phase point passes another with a finite difference in their velocity. Because 
the relative velocity has no singularity as one particle passes the other, the PDF of the 
spatial separation also has no singularity. It may therefore be approximated by a uniform 
distribution in the vicinity of dx = 0. The corresponding distribution for Y is obtained 
by writing the probability element as follows dP = ps x (dx)ddx = py(Y)<1Y. Hence 

Py{Y) ~ const, x -jpr ~ exp(T) (20) 

This contribution is negligible compared to that from the analysis of the differential equa¬ 
tion whenever the latter predicts a < 1. The contribution from the folding events is 
therefore smaller than that due to fractal clustering whenever D 2 < 1. This condition 
is never violated in one dimension | 2 T]. In higher dimensions the equation analogous to 
(1151) does not have finite-time singularities, although there are caustic singularities where 
volume elements vanish 
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3 Markovian approximations 


We wish to use information about statistics of Z\{t) to determine Z ?2 = ot- The most 
practicable approach is to use a Markovian assumption, where the future development of 
a system can be assumed to be independent of its past history. In the present context, 
we assume that future evolution of Y{t) is determined by its current value, and by the 
current value of Z\. We therefore consider a joint PDF of Y and of Z\. Given Y and Z i, 
let K(AY, Z\, Z[,t ) be the PDF for Y to increment by Ay and for Z\ to reach Z' x after 
time t. The joint PDF of Y and Z\ evolves according to 

/ OO pOO 

dAY / dZ[ K(AY, Z[,Z X , At) p Y , Zl (Y — AY, Z[,t — At) . (21) 

-oo J — OO 

The steady-state probability density is expected to be a product 

PY,z i (Y, Z x ) = p Zl (Zi) exp(ay) (22) 

where the distribution p Zl (Z\) will be discussed shortly. Because equation (12T1) is derived 
by linearisation of the equations of motion, equation (I22j) is valid in the limit as Y —>■ —oo. 
In order for the distribution to be normalisable, we require that p Y , Zl approaches zero 
sufficiently rapidly as Y —>• —oo, implying that a > 0. 

Inserting (1221) into f[2lj) . the steady-state distribution p Zl (Z i) and the exponent a must 
satisfy an integral equation 

/ oo poo 

dAY dZ[K{AY,Z[,Z l ,At)p z pz[)eM-^Y) (23) 

-OO J —OO 

which is valid for all At. 

Consider the distribution p Zl (Z\) in (1221) . It might be expected that this is the same 
as the distribution of Z\(t) obtained from equation (|18jl or its multi-dimensional general¬ 
isation. We term this distribution po(Z\). However, the distribution p Zl (Z\) differs from 
Po(Zi) because it is conditioned upon being at a particular value of Y [20]. If a / 0, 
particles reaching a negative value of Z\ have recently arrived from a larger value of Y, 
where the probability density is larger. This implies that the distributions are different, 
and moreover that the distribution p Zl (Z\) has a smaller mean value than po(Zi). 

Now consider the application of this equation in two limiting cases. 

3.1 Short propagation time 

Consider the limit At —>• 0 in (|23|h In this limit the structure of the propagator can be 
simplified, because AY = ZAt + 0(At 2 ). This implies that one of the integrals can be 
eliminated from (1231) . and we may write 

/ OO 

dZ[ W(Z{,Zi, At)p Zl (Z0exp(-aZ{At) (24) 

-oo 

where U(Z[, Z\, At) is the propagator for the random process Z\[t) with equation of 
motion m (or its higher-dimensional generalisation) to reach Z\ from in time At. 

If a Markovian approximation is valid in the limit At —>• 0, we have continuous-time 
Markov process for ([23l) . and the probability density p Zl (Z\,t) obeys a Fokker-Planck 
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equation [25], where the evolution kernel U(Z[, Zi,t) is generated by a Fokker-Planck 
operator F: 

dj ^~ = FPZ, ■ (25) 

We can represent functions as vectors using Dirac notation, so that (f25l) is notated as 
follows: 

dtlpzj = P\pz 1 ) ■ (26) 

For small values of At the action of the propagator kernel can then be approximated by 
Li (A t) = Z + F At + 0(At 2 ), where 1 is an identity operator, that is for a function f(Z\) 
represented by a vector |/), we have 

/ OO 

dZ[U(Z[,Z 1 ,At)f(Z[) = \f)+F\f)At + 0(At 2 ) . (27) 

-oo 

For small values of At equation (1241) then reduces to 

pzAZi) = exp(-aZ 1 At)p Zl (Zi) + AtFp Zl {Z\) + 0(At 2 ) . (28) 

In the limit as At —>• 0, this relation implies the condition 

[F — aZi] pz 1 (Z 1 ) = 0 (29) 

which is a partial differential equation for pz 1 (Z \) and a. 

At this stage it is useful to consider a concrete example. The one-dimensional model 
equation of motion for Z \, Eq. (fl8l) . can be regarded as a stochastic differential equation, 
in which the velocity gradient E(t) is a random element. If the correlation time of E(t) is 
sufficiently small, a Markovian approximation is justified, and E{t) can be replaced by a 
multiple of a white noise signal, p{t) 1 which has the following statistical properties: 

(v(t)) = 0 , {r]{t)v(t')) = 5(t - t') . (30) 


The equation of motion for Z\ is replaced by 

Z\ = — 7 .Z 1 — Z\ + VzDpit) (31) 

where the diffusion coefficient is 

V = Y fZ dt{mEm ' (32) 


The Fokker-Planck operator corresponding to the Langevin equation (1311) is |29j : 

/=( 7 Z 1 + Z?)A_ +p J4 (33 ) 

so that for the one-dimensional model equation (I29[) reduces to an ordinary differential 
equation 


d 

dzl 


{ 1 Z l + Zl)pz 1 (Zi) + ^(Z 1 ) 


ol Z 1 pz 1 (Z 1 ) = 0 . 


(34) 
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We require normalisable solutions pz t (Zi), which only exist for particular values of a. 
(Later, we give a prescription leading to a unique series solution of this equation). Upon 
integrating over space, and using the fact that T is a divergence, we have 



dZ\ Z\ pz 1 (Zi) 


(Zi) = 0 


(35) 


The equations (12171) and (1351) together constitute a new and exact method for determining 
L >2 = a, in two steps. First equation (1291) is solved to determine a one-parameter family 
of solutions. Second, the correct value of Z ?2 is determined by finding the value of a for 
which the mean value of Z\ is zero [20]. 

This approach has the attractive feature that it involves the analysis of differential 
equations, which are susceptible to many types of mathematical techniques. 


3.2 Long-time propagation 

In the long-time limit we expect that a Markovian approximation is always valid. In this 
limit the propagator is expected to ‘forget’ the initial distribution, so that 

K(AY, Z[, Z l5 At) = p Zl (Z[) pay (AY, At) (36) 

independent of Z\ , where pay (AY, At) is the probability of a displacement Ay in time 
At. 

We now apply the large-deviation principle EE EU to the statistics of AY. This 
principle concerns the statistics of time averages such as the finite-time Lyapunov exponent 
a(t) = AY/t, equation (fTP) . It is expected that the tails of the distribution pay (AY) 
satisfy 

Pay (AY, t) ~ exp[-tI(AY/t)} (37) 

for some function 1(a), which is termed the ‘entropy function’ in the literature on large- 
deviation theory. 

The displacement is Ay = at. Changing the variable of integration in (123[) from Ay 
to a, we obtain 

/ oo roo 

da dZ[ pz 1 (Z[)p Zl (Zi)ex-p[-t(I(a) + aa)\ . (38) 

-OO J — OO 

Assuming that pz 1 (Z\) is a normalised distribution, this gives 

/ OO 

d <7 exp[— t(I(a) + aa)} . (39) 

-OO 

This integral is an implicit relation between a and the large deviation function 1(a). The 
integral is estimated using the Laplace principle: in the limit as t —>• oo, the integral is 
estimated by determining the value of the integrand at its maximum. The maximum is at 
position a* determined by the condition 

I'(a*) = -a (40) 
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and the integral is estimated as t exp[—t(I(a*) + ua*)\ ~ 1, so that 

I{a*) + aa* = 0. (41) 

These equations, (1401) and (1411) can, in principle, be solved to determine a = D 2 . Similar 
approaches are discussed in [is m Ei- The difficulty lies in determining the entropy 
function, I(cr). 


3.3 An approximate expression for D 2 

Before exploring the applications of equation (l29l) in greater depth, we describe an ap¬ 
proximate expression for D 2 , previously discussed in [32], which is asymptotically correct 
in the limit as D 2 —>• 0. Because of its simplicity, it is a natural benchmark against which 
other approaches can be compared. 

We observe that the variable Y has diffusive fluctuations 


((A Y(t) - At) 2 ) = 2VyAt 


(42) 


with diffusion coefficient 

/ OO 

dt [(Z 1 (t)Z 1 (0))-(Z 1 ) 2 ] . (43) 

-OO 


On time scales which are large compared to the correlation time of Z\(t) we expect that 
the probability density P(Y, t ) satisfies a Fokker-Planck equation 

dP d <9 2 

¥ = -w ( ' p ) + ^ p »' (44) 

The drift velocity v and diffusion coefficient T>y are defined by the relations 


v = 


m 

St 


Vy 


(<^ 2 ) 

2 St 


When v and T>y are constant, this equation has an exponential solution 


(45) 


P = Aex.p(^-Y^J . (46) 

Noting that the drift velocity v is equal to the Lyapunov exponent, (Z i) = A, comparison 
with (fTOl) and (fl2l) implies that 

D 2 = a = A . (47) 

This approximation is only valid when A > 0, because no normalisable solution can be 
constructed if P(Y ) is diverging as L —>■ — oo. 

The use of the Fokker-Planck equation is only justified when the gradient of P(Y,t) 
is sufficiently small. The condition is that dP/dY should be small compared to \/8Yq, 
where 8Y$ is the scale over which Y varies during its correlation time. The condition for 
the validity of (1471) is therefore (Z\)/(\Z\\) <C 1, which is equivalent to D 2 <C 1. 

Consider how equation (|47l) relates to the long-time limit of the propagator. The statis¬ 
tics of the displacement Y(t) are directly related to the finite-time Lyapunov exponent: 
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A Y(t) = a(t). The variance of a(t) is 2 Vyt. In the case where 1(a) can be adequately 
approximated by a quadratic function, we see that 1(a) may be approximated by 


1(a) 


(o- - A ) 2 

4 Pi- 


Using this approximation in (1401) and (1411) we recover Eq. (1471) . 


(48) 


4 Three-dimensional model 

In this section we consider how to compute the correlation dimension for inertial particles 
suspended in a three-dimensional flow. In order to make it possible to perform the analy¬ 
sis, we consider particles in a random velocity field with known statistical properties. This 
approach has been successful in modelling the Lyapunov exponents of particles in turbu¬ 
lent flows: the leading Lyapunov exponent was obtained in [33], and all three Lyapunov 
exponents for the spatial separation of particles in [25], showing excellent agreement with 
the numerical simulations of particles in turbulent flows described by Bee [7]. Here we 
build upon the results of these earlier calculations by analysing the correlation dimension 
for the same random-flow model. 

The flow underlying turbulent aerosols is usually incompressible, V • u = 0. But in 
order to analyse the properties of the perturbation theory employed in this paper, it is of 
interest to also consider partially compressible flows. We use the following decomposition 
of the flow velocity into solenoidal and potential contributions 

M = C 3 (VAA + fiVip) (49) 

where C 3 is a constant. This model was used in [33] to compute the maximal Lyapunov 
exponent of inertial particles in random, partially compressible flows. The parameter 
/3 determines the relative magnitude of the potential and solenoidal contributions. A 
convenient measure of the relative importance of these two contributions is 

r _ 4 + / 3 2 
2 + 3/3 2 

Since the parameter /3 assumes values between zero and infinity, we have that ^ < T < 2. 
The case T = 2 corresponds to solenoidal flow (/3 = 0). For T = |, by contrast, the flow 
is purely potential ((5 —>• 00 ). A special case of interest discussed below corresponds to 
T = 1, where the solenoidal and potential contributions are of equal strengths. 

We take the components of A and ip = Aq to be Gaussian homogeneous isotropic 
random functions with zero mean values and correlation functions 

(A i (r,t)A j (r',t')) = 8ijC(\r-r'\,\t-t'\) . (51) 

The correlation function C is assumed to decay to zero for spatial separations much larger 
than the correlation length 77 of the flow, and for time differences much larger than the 
correlation time r. The typical fluctuation size of the flow is denoted by ( u 2 ) = v,q. This 
implies that the normalisation constant in (|49|) must be chosen as: 

Ci = u 2 o [3(2 + p 2 )\C"(0,0)\]~ 1 . (52) 



11 




Following the approach in [33] and [35] we analyse this model in the ‘white noise’ limit 
r —» 0, which justifies the use of the Markovian approximation considered in section 13.11 
The fluctuations of the velocity gradients E(t) are characterised by specifying a set of 
diffusion coefficients, analogous to equation (1321) . The diffusion coefficients are expressed 
in terms of the correlation functions of the elements of E: 


Du = 


t r 
2 


df (E a (t)En(0)) 


(53) 


(the factor of y 2 is a consequence of the fact that E(t) is multiplied by 7 in (1TH1) ). There are 
some technical complications involved in calculating the Fokker-Planck operator appearing 
in (1291) . which were discussed in detail in [33]. We can read off the Fokker-Planck operator 
from the results in that paper. The version of (1291) which is applicable to our model is 


dp 

dt 


= Jzi^ 21 + z > “ z * - zi)p\+ 

+ 5^ [(7Z2 + 2ZlZ2)p]+P22 l4 

d a 2 p 

+ Tj7r[(7-^3 + 2ZiZ 3 )p] +T>33^2 ~ aZip. 


(54) 


The steady-state form of this equation is analogous to the one-dimensional equation (|33j) . 
In three dimensions the condition (|35l) takes the form 

/ OO POO POO 

dz 1 / dZ 2 / dZ 3 Zi p(Z ll Z 2 , Z 3 ) = (Zi) = 0 . (55) 

-00 J— OO J— 00 

The correlation dimension (equal to a) is obtained by finding a value of a for which a 
normalisable solution of (|54D can be obtained for which the mean value of Z\ is zero. 
The equations (154[) and (1551) together constitute an exact method for determining the 
correlation dimension in the white-noise limit. 


5 Perturbation theory 

Here we derive a perturbation expansion for the correlation dimension. It is convenient to 
introduce dimensionless variables: 


Xi = \Jl/Da Zi. (56) 

The expansion parameter of the perturbation expansion is given by e where 

2 _ ^11 
e “ 7 3 

Because e —» 0 in the over-damped limit, this perturbation parameter plays a role which 
is analogous to the Stokes number. The connection between e and St is discussed in 
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detail in [2S]• We denote the joint probability density of x 3 ,..., x 3 in the steady state by 
P(x i,x 2 ,x 3 ). It follows from Eq. (1541) that P satisfies the equation: 


0 = PP=-^-[( Xl +e{xl- r(x| + x|)))P] 

Q Q 

+ ^-[x 2 (l + 2exi)P] + — [x 3 (l + 2exi)P] 

d 2 P d 2 P d 2 P 

+ o 2 o 9 + o 2 _ 6ax iP • 

oxf ox 2 axg 


(58) 


This equation defines the Fokker-Planck operator F(e, a, T). Following {2(TI . we now de¬ 
velop its solution as a series expansion in e, using a system of annihilation and creation 
operators which are analogous to those used in quantum mechanics. We employ a nota¬ 
tion similar to the Dirac notation: a function /(xi,X 2 ,X 3 ) is denoted by a vector |/). The 
scalar product between two states |/) and | g) is given by 

poo poo poo 

(f\g)= dxi / dx 2 / dx 3 f{x I,x 2 ,x 3 ) g{x i,x 2 ,x 3 ). (59) 

J— OO J —OO J— OO 


We expand both the solution |P) of (|58l) and the value of a for which the solution of this 
equation exists and satisfies (xi) = 0 as power series in e: 


fc =0 fc =0 

The Fokker-Planck operator in Eq. (|58D is written as 

F = Fq + e(C/ - axi) 


(60) 


(61) 


which defines the operators Fq 


and Q. The unperturbed steady-state |Pq) satisfies 


Fol-Po) = 0. 


(62) 


It is given by 


Pq{xi,X2-,xP) 


exp[— {x\ + xj + x§)/2] 
(27t) 3 / 2 


(63) 


Other eigenfunctions of Fo are generated by creation operators di and annihilation oper¬ 
ators bi: 


di — dxi 

bi = d Xi + Xi . (64) 

These operators generate eigenfunctions satisfying 

F"o| (frpnm) — (jl + Til, T p')\(f > pnm') (65) 
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according to the rules 


^1 |0p,n 

,m) 

— \ c t > p+l,n,m) 


m) 

— P\&p—l,n,m) 


,m) 

= \$p,n+l,m) 

^11 

m) 

— ^\^p,n— l,m) 

^31 0p,7l 

,m) 

= \4 > p,n,m+l) 

^11 

m) 

= m\(j)p jntrn -i) 


(66) 


with |<^ooo) = \Po), normalised as a probability density. The states | Pk) in (1601) are 
expressed as linear combinations of the eigenfunctions 1</> pnm ): 


oo oo oo 

P k) = S S S Ppnm l&mm) • (67) 

p=0 77/—0 ra=0 


The eigenfunctions generated by repeated applications of a are neither normalised nor do 
they form an orthogonal set. This is different from earlier perturbation theories for the 
Lyapunov exponent of inertial particles [HU !25] . 

We first consider how the condition (aq) = 0 constrains the coefficients Ppnm in ([BTJh 
Using (IMP we find: 

(®1> = (zil-P) = e * Ppnm 11 4*pnm) = > (68) 

k p,n,m k 


so that the condition (aq) = 0 is satisfied by requiring that 

Pi oo = 0 (69) 

for all values of k. Substituting (1601) into (|61li leads to a recursion for | Pk). The term of 
order e k is given by 


fc-i 


0 = Jo I Pk) + G\Pk-l) — ^ oq(ai + &i)|-Pfc-j_i), 

1=0 


(70) 


with 


and 


Jo = —a\b\ — a 2 6 2 — ^3^3 


(71) 


Q = -ai[((ai + h ) 2 - T((a 2 + b 2 f + (a 3 + & 3 ) 2 ))] 

— 2(ai + &i)[a 2 (a 2 + S 2 ) + 03(03 + £>3)] • (72) 

Equation (1701) is a recursion for the state | Pk), and the coefficient ak -1 in terms of the states 
| Pj) and coefficients aj determined in previous iterations. By considering the coefficient 

of 1 4> P nm), we obtain Ppnm in terms of the coefficients p'^, J\ n , and oq, with l = 0,. .., k — 1 . 
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In order to determine the coefficients cq, consider the case p = n = m = 0, where the 
coefficient of the state |0ooo) i n Eq. (1701) reduces to the condition: 

S>*&' -1) = 0- (73) 

1=0 

This condition can be fuffilled in at least two ways. One solution is obtained by setting 
all ctk equal to zero. This case corresponds to calculating the Lyapunov exponent. Using 
A = (Z\) = sJVn/'y (x\) we find Eqs. (67,68) in [33]. A second possibility is to require 
(xi) = 0 corresponding to p) 0 g = 0 for all values of k, as explained above. This is the case 
relevant for calculating the correlation dimension. Using the initial condition 

Ppnm = ^pO^nO^mO (74) 

(fc) 

we can iterate Eq. (I7U|) to determine the coefficients Pp nm hr terms of the ai . The first few 
non-vanishing coefficients obtained by recursion of (17011 are listed in Table CD Note that 
Eq. (1701) only relates coefficients Ppn m with indices n and m to coefficients with indices 
n 1 and m! provided that v! — n and m! — m are even integers. This implies that only 
coefficients with even values of n and m are non-zero (see Table [T]). Note also that if all 
odd-order coefficients 0 ( 2 n+i vanish, then Eq. (|70l) does not mix the parity of p + k in ppnm■ 
We find that in this case, Eq. (1701) provides two independent recursions: one for ppn m with 
even p + k and the initial condition p 2 p n m = ^po^no^mo when k = 0, and one for ppnm 
with odd p + k and the initial condition P 2 p+i n m = 0 when k = 0. Using the boundary 

condition p f^ = 0 (Eq. ([75l) l we find that indeed all odd order a 2 n +i vanish, and that all 
coefficients with odd values of p + k vanish. This is illustrated to the lowest order in k in 
Table CD all displayed coefficients which are of different parity in p and k are multiplied 
by a\ = — J>ioo = 0. 

These considerations do not determine the normalisation of the distribution |P). Ex¬ 
panding the normalisation condition in terms of Eqs. (|60l) and (|67l) yields 

Pooo = 4o • (75) 

Once the coefficients Pmhp have been determined to each order k. we use Eq. (1691) to 
compute a.k- 1 - From Table CD we find the first three coefficients in the expansion ()60D of a 
in powers of e: 


a 0 = 2T —1, «i = 0, a 2 = -2(r-l)T(2r + l). ( 76 ) 

This gives the correlation dimension for the three-dimensional random-flow model in the 
white-noise limit 

D 2 = 2T - 1 - 2T(T - 1)(2T + l)e 2 + ... ( 77 ) 

to second order in e. As mentioned in the Introduction, the two leading non-zero coef¬ 
ficients of this expansion for incompressible flows (T = 2) were computed in [19] . The 
coefficient of e 2 in that work differs from our result. 
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6 Results and discussion 


We use an algebraic manipulation program to obtain the series expansion of D 2 (e) in 
powers of e from Eq. (1701) to higher orders in k. To order e 10 the result is 

D 2 = 2T - 1 + 2T(T - 1)(2T + 1)[—e 2 + (-11 -2T + 6T 2 )e 4 
+ ± (-588 + 5T + 391T 2 - 16T 3 - 144T 4 )e 6 

+ ±(-42579 + 18573r + 22727T 2 - 19284P 3 - 12648T 4 + 3464T 5 + 3960T 6 )e 8 
+ ^ (-3863052 + 3918303T + 288351T 2 - 4153120T 3 + 47186T 4 

+ 1409736T 5 + 277928T 6 - 216448T 7 - 117936r 8 )e 10 ] 
+ 0(e 12 ). (78) 

In special cases we have obtained expansions to higher orders. For incompressible flows 
(T = 2), for example, we find to order e 28 : 

£>2 = 3- 20e 2 + 180e 4 - 9640e 6 + 206940e 8 

- 16548920c 10 + 477315000c 12 - 50149424368c 14 
+ 1692947357004c 16 - 5614110582647928/25e 18 
+ 209543657412608424/25e 20 

- 860424252594210743568/625e 22 

+ 241528608428504721258888/4375e 24 

- 8471768050800513607578954992/765625e 26 

+ 2514450499347358305045823304592/5359375c 28 
+ ... (79) 


For incompressible flows we have obtained the first 33 non-vanishing coefficients as frac¬ 
tions of integers and the following 17 coefficients to ten significant digits. 

The corresponding series in two spatial dimensions was derived by [20] : 

£>2 = T - 1 - T(T 2 - l)e 2 + r(T 2 - 1)(3T 2 + 2 T - ll)e 4 + 0(e 6 ). (80) 


Iterating the recursions derived by m we have obtained the first non-vanishing 110 coef¬ 
ficients to ten significant digits in the incompressible case (T = 3). 

The series quoted above are asymptotically divergent: they diverge but every partial 
sum of the series approaches D 2 as e —)• 0. Evaluating the coefficients a k for a given value 
of T shows that they grow factorially as a function of k: 

a k ~ a r S^ k (k - 1)!(1 - b v /k + •••)• (81) 


This is a typical asymptotic behaviour of the coefficients a k for large values of k (35] . 
Here the ‘action’ Sr and the constant br are obtained by fitting of the ansatz (1811) to 
the coefficients. For the fit we use a non-linear least-squares method, assuming that the 
relative magnitude of subsequent coefficients is on the form 


a k+ i ! k 2 k + 1 - br 

- rs -' Op- 

a k 1 k + 1 k — br 


(82) 
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Figure 1: ( Color online), a Fit to asymptotic form (18211 (red solid line) for T = 3/2 and 
d = 3 in the range 25 < k < 50. The quotients —a^+i/ctk are shown as blue symbols. 
The resulting values of the fitted parameters are Sr=3/2 = 0-167 ^r= 3/2 = 4.8. b Shows 
actions obtained from the perturbation coefficients as a function of T, in two and in three 
spatial dimensions. The first 110 non-vanishing coefficients in two dimensions and the 
first 50 non-vanishing coefficients in three spatial dimensions are fitted to (1821) using a 
non-linear least-squares method. The fitted actions Sr are shown as symbols: red o in two 
dimensions, and green □ in three dimensions. Also shown are the curves S' = 1/6 (dotted 
black) and S = 1 /(6|T — 1|) (dashed black). The coloured regions indicate the allowed 
ranges of the parameter T: 1/3 < T < 3 in two spatial dimensions (red) and 1/3 < T <2 
in three spatial dimensions (hashed green). 


Fig. dk illustrates the asymptotic behaviour of the coefficients, using the case T = 3/2 in 
three spatial dimensions as an example. The action Sr extracted from fits such as the one 
in Fig. [Br is shown in Fig. db^ in both two and three spatial dimensions. The resulting 
action is found to depend upon T as follows: 

S r = min[l/6, l/(6|r - 1|)], (83) 

in both two and three spatial dimensions. We note that the coefficients of the perturbation 
series for the maximal Lyapunov exponent in two spatial dimensions [3j] give rise to the 
action 1/(6|T — 1|) for all values of T. 

We also note that the two- and three-dimensional cases shown in Fig. db differ from 
each other. In three dimensions the action is always given by 1/6 in the allowed range 
of r (this is not the case in two spatial dimensions). As opposed to the perturbation 
expansions for the Lyapunov exponent and the two-dimensional correlation dimension, 
the three-dimensional perturbation expansion for the correlation dimension is determined 
by one action only, S = 1/6. 

We have resummed the perturbation series (17811 and (18011 using Pade-Borel resumma¬ 
tion: to sum the series 

OO 

T>2(e 2 ) ~ a 2/e 2 *, (84) 

l=o 

consider the modified series, the so-called ‘Borel sum’ (assumed to have a finite radius of 
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convergence due to the extra factor of 1 //!) 

oo 

^ 2 ) = E ^ w 

1=0 L ' 

Then the sum is estimated by 

D 2 (e 2 ) = Re f d te~ t B(e 2 t). (86) 

Jc 

The integration path C is taken to be a ray in the upper right quadrant of the complex 
plane. In order to perform the integral, an approximation of the Borel sum outside its 
radius of convergence is required. One possibility is to approximate B by ‘Pade approx- 
imants’ [36] of order [n,n] (or [n,n+ 1]). For T = 2, the Pade approximations of order 
[2,2] and [3,3] are (with x = e 2 ): 

, , x 2 -20x , 

-^[2,2] ( x ) — 3 + 1671 649927 2 ’ (^7) 

1 1442 x 8652 x 

728234642879 „3 , 91933567500 on 

R / \ _ o i_ 562607103 _ ~r 187535701 x x 

^[3,3] W o ~r _ 7505535441 _ 99077893373 2 , 11726142610857 3 ' 

1 375071402 x 937678505 x _+ " 7501428040 x 

Higher orders are too lengthy to write down here. Fig. [2] shows the results we obtained for 
D 2 for T = 2 in three spatial dimensions by integrating H[ 4 4 ], _B[ 8) 8], and H[ 1616 ] according 
to Ecp (1861) . The corresponding contour C in the complex f-plane was chosen along a ray 
from the origin at angle 7r/4. For small values of e the results depend only negligibly on the 
precise choice of the contour. Also shown are results for Z ?2 obtained by direct numerical 
simulations of the equations of motion Q . We observe that the Pade-Borel resummations 
converge quickly for not too large values of e, and we find excellent agreement with results 
of direct numerical simulations of the random-flow model. 

It is clear, on the other hand, that the resummation fails for larger values of e. We 
suspect that a non-analytical contribution of the form A exp[—l/(6e 2 )] is not captured and 
must be added to the perturbation series. In [34] it is shown that a corresponding term 
must to be added to the perturbation result for the maximal Lyapunov exponent. The 
situation here is similar. This is most easily seen by considering the case T = 1. Eq. (1781) 
shows that the first twenty perturbation coefficients vanish for T = 1. We hypothesise that 
all coefficients vanish at T = 1 and that the correlation dimension exhibits a non-analytic 
dependence on e, of the form 

D 2 ~ A 1 exp ( - . (88) 

This is shown in two and three spatial dimensions in Fig. [3] These results complement 
earlier studies of the information dimension D 1 , discussed in detail in [25]: the Borel 
summation technique was more successful in that case, but D\ has less direct physical 
significance. 

We conclude this section with two further comments. First, for small values of e we see 
that in incompressible flows 3 — D 2 oc e 2 oc St. This is a consequence of the fact that we 
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Figure 2: (Color online). Correlation dimension for the white-noise model in three spatial 
dimensions for T = 2 as a function of e. Shown are results of direct numerical simulations 
of the equation of motion ()16l) . symbols, and results of Pade-Borel resummations of the 
perturbations series for D 2 , of order [4,4] (dash-dotted line) [8,8] (dashed line), [16,16] 
(solid line). 


considered the white-noise limit. In this limit the fractal information dimension exhibits 
the same scaling [25]. In flows with finite correlation time, by contrast, the correlation 
dimension deficit behaves as 3 — D 2 oc St 2 for small Stokes numbers [371 EHl |39]. 

Second, we note that the correlation dimension exhibits a singularity in the advective 
limit (e = 0) as the compressibility parameter approaches T = 1/2, corresponding to a 
path-coalescence transition where the maximal Lyapunov exponent changes sign [33]. The 
perturbation theory gives correct results for T > 1/2, it fails for T < 1/2. 

7 Conclusions 

In this paper we have derived a general method for calculating the correlation dimension of 
random dynamical systems, which complements DNS (direct numerical simulation) studies 
of particles in turbulence BIS] and numerical studies of stochastic models [22, [23], [25] . The 
method is formulated in terms of a propagator describing the time evolution of particle 
separations and particle-velocity gradients. In special cases, known methods for computing 
the correlation dimension are obtained [Ml H6l [31132]. 

A short-time expansion of the propagator yields a solvability condition on a partial 
differential equation, leading to a perturbative expansion of the correlation dimension, for 
which the coefficients can be obtained exactly and to any order. We derived the exact 
first 33 coefficients in a series expansion of the correlation dimension for inertial particles 
in three-dimensional spatially smooth random flows that are white noise in time. Related 
series expansions have been presented for Lyapunov exponents of inertial particles in such 
flows in earlier works [33] [251 El] 

We have obtained accurate results for the correlation dimension of inertial particles in 
three-dimensional white-noise flows by Pade-Borel resummation of the perturbation series 
for not too large values of e. However, for the correlation dimension D 2 the resummation 
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Figure 3: ( Color online). Correlation dimension for T = 1 as a function of e~ 2 in two spatial 
dimensions (a) and in three spatial dimensions (b). Also shown is the non-analytical law 
(|88l) with prefactors A\ = 1 in two dimensions and A\ = 2 in three dimensions, (dashed 
lines). 


method is not as successful as for the information dimension D\, which was considered 
in [25]. It would be desirable to develop a more direct analytical approach to extracting 
information about D 2 from equation (|54lh 

In a particular case, for T = 1, we find that the perturbation coefficients vanish and 
the correlation dimension exhibits a non-analytical dependence upon e. We conjecture 
that there is a corresponding non-analytical contribution also for T > 1. 

Finally, we remark that it is possible to extend the method presented here to treat ve¬ 
locity fields with finite correlation time. This can be achieved by considering a temporally 
smooth velocity gradient obtained from a stochastic process which is driven by a white 
noise signal. Numerical studies of velocity gradient statistics in turbulence show that they 
have correlation functions which are well approximated by exponentials [40]. Velocity 
gradients of turbulent flows can, therefore, be modelled by an Ornstein-Uhlenbeck process 
m, as described in m- The operator methods used in this present paper have been 
extended to temporally smooth velocity gradients m, but some care may be required in 
their application and interpretation [43] [43] . 
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Table 1: Non-zero expansion coefficients in Eq. fl67l) for k = 0,..., 3. Obtained by 
recursive solution of Eq. <©■ 


(o) _ i 

P 000 — 1 


Pioo — 2 T 1 ao , P102 — P120 — F 3 2 > P300 — 3 


„(2) _ (2) _8T^7 , r^8 (2) _ £^2 (2) _ (2) _ 1^2 (2) _ (4T-5)(4r~3) . 4r~5 2 /g 

”002 — 1*020 — 3 ^ 6 a ° ’ -C022 — 3 ’ *3)04 — 1*040 — 6 ’ ^200 — 6 ^ 2 “0/ z > 

„(2) _ ( 2 ) _ (-13+17r-6r 2 ) r-2 (2) _ ( 2 ) _ ffi- 2) 2 ( 2 ) _ (r- 2) 2 ( 2 ) _ 4T^5 _ 20 „0) _ _(2) _ r ^2 

"202 "220 6 I 3 ^0 ? ^204 "240 18 ’ **222 9 ’ "400 6 3 ’ ^402 "420 9 ’ 


r 


pfdo = -5 + 20T - 16 r 2 + 2(— 5 + lor - 3r 2 )a 0 + 2(2r - 3)a§ - ag - a 2 , 

„(3)_(3) _ -130+247r-118r 2 , -20+19r-2r 2 , 5r-28„,2 „(3)_(3) _ -122+167r-56r 2 , -22+13r-r 2 

1*102 — P 120 ~ 9 "T" 2 a 0 -I- 18 «o , p 104 — P 140 — 30 ~r 18 


(3) _ -122+167r~56r z _|_ -22+13r~r 2 


P 122 — 


15 


(3) _ (3) _ (r- 2 ) 2 


-ao , p 124 — P 142 — 


„(3) _ -135+290r-212r 2 +56r 3 , -67+77r-24r 2 , 9r-19„,2 

P 300 “ 18 1 9 a ° ^ 9 “ n 

„(3) _ „(3) _ -715+1194r-724r 2 +172r 3 , -65+6ir~18r 2 
1*302 — 1*320 — 90 L8 


6 ’ 
3 

iffl 


0 6 ’ 

r-2 2 J3) _ J3) 

~E~ a 0 ; 


1*304 = 1*340 = §(—6 + 9T - sr 2 + r 3 ) - 
P% 2 = |(-6 + 9T - 5r 2 + r 3 ) - l/9(—2 + r) 2 a 0 , pfl = pS 2 = , 


(r- 2 ) 2 
18 a ° 































